home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / bit / src / holo.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  10KB  |  437 lines
  1. /*
  2.  * $Id: holo.c,v 0.91 1994/02/20 00:52:53 zhao Pre-Release $
  3.  *
  4.  *. This file is part of BIT shareware package. After the two weeks of
  5.  *  free evaluation period, you are encouraged (required) to register
  6.  *  your copy for a small registration fee, which is $35 for personal use
  7.  *  and $50 for commercial, government and institutional use.
  8.  *
  9.  *  Copyright(c) 1993, 1994 by T.C. Zhao.
  10.  *  All rights reserved.
  11.  *
  12.  *  Permission to use, copy, and distribute this software in its entirety
  13.  *  for non-commercial purposes is hereby granted, provided that the
  14.  *  above shareware and copyright notices and this permission notice
  15.  *  appear in all copies and their documentation.
  16.  *
  17.  *  This software may be modified for your own use, but modified versions
  18.  *  may not be distributed without prior consent of the author.
  19.  *
  20.  *  This software is provided "as is" without expressed or implied
  21.  *  warranty of any kind.
  22.  *
  23.  *.
  24.  *
  25.  * Read in X-ray Photo-emission data file and generate intensity maps
  26.  * (Hologram)
  27.  * default: Raster is assumed to be from LEFT to RIGHT, TOP to
  28.  * BOTTOM. Aspect ratio will be kept as the input.
  29.  *
  30.  * To force Xsize==Ysize, use @ASP. To force top to bottom and LEFT to
  31.  * right(rot90), use @YX To force Bottom to TOP , use @B2T (rot180)
  32.  *
  33.  * Note:
  34.  *
  35.  * All matrix are arranged as Matrix[Y][X], i.e.; X runs first
  36.  */
  37. #if !defined(lint) && defined(F_ID)
  38. char *id_holo = "$Id: holo.c,v 0.91 1994/02/20 00:52:53 zhao Pre-Release $";
  39. #endif
  40.  
  41. #include "bit.h"
  42.  
  43. #ifdef ADD_HOLO
  44.  
  45. #include "extern.h"
  46. #include "dmalloc.h"
  47.  
  48. static void holo_getmap(IPTR, FILE *, int);
  49. extern void *interpol2(ci_t **, int, int, int, int, int);
  50. static int raw, cmin, cmax;
  51. static char tmpstr[120];
  52. static xyequal = 0, yx = 0, t2b = 1, mag = 0;
  53.  
  54. static int lindex = 1;
  55.  
  56. int
  57. HOLO_desc(IPTR im)
  58. {
  59.     int i, j, err = 0;
  60.     int xs, ys;
  61.     FILE *fp = im->fp;
  62.  
  63.     if (Badfgets(tmpstr, sizeof(tmpstr) - 1, fp))
  64.     return -1;
  65.     raw = (strncmp(tmpstr, "DSBB", 4) == 0);
  66.  
  67.     /* dimension must follow */
  68.     fgets(tmpstr, sizeof(tmpstr), fp);
  69.     if (sscanf(tmpstr, " %d %d ", &xs, &ys) != 2)
  70.       {
  71.       Bark("HOLO_desc", "Error getting desc");
  72.       return -1;
  73.       }
  74.     while (fgets(tmpstr, sizeof(tmpstr), fp) && tmpstr[0] == '#')
  75.     ;
  76.  
  77.     err += (sscanf(tmpstr, "%d %d", &i, &j) != 2);
  78.  
  79.     while (fgets(tmpstr, sizeof(tmpstr), fp) && tmpstr[0] == '#')
  80.     ;
  81.     err += (sscanf(tmpstr, "%d %d", &i, &j) != 2);
  82.  
  83.     while (fgets(tmpstr, sizeof(tmpstr), fp) && tmpstr[0] == '#')
  84.     ;
  85.     err += (sscanf(tmpstr, "%d %d", &cmin, &cmax) != 2);
  86.  
  87.     while ((i = getc(fp)) != EOF && i == '#')
  88.       {
  89.       while (getc(fp) != '\n')
  90.           ;
  91.       }
  92.     ungetc(i, fp);
  93.  
  94.     im->w = xs;
  95.     im->h = ys;
  96.     im->colors = im->cmap->colors = (cmax - cmin + 1);
  97.     return err ? -1 : 0;
  98. }
  99.  
  100. /* for input size smaller than {hw}_min, interpolate to def_{hw} */
  101. static int w_min = 201, h_min = 201;
  102. static int def_w = 221, def_h = 221;
  103.  
  104. #include <stdlib.h>
  105. int
  106. HOLO_load(IPTR im)
  107. {
  108.     int i, j, err = 0;
  109.     float fcc;
  110.     int nh, nw;
  111.     int xfac, yfac, factor;
  112.     ci_t **pp, **pp1;
  113.     FILE *fp = im->fp;
  114.     float xunit = 1.0, yunit = 1.0;
  115.     int xori = 0, yori = 0;
  116.     long rlines;
  117.  
  118.     rlines = progress_report("Loading Hologram ...", im->h);
  119.     pp1 = im->mraster;
  120.     if (raw)
  121.       {
  122.       for (j = 0; j < im->h; j++)
  123.         {
  124.         REPORT(j, rlines);
  125.         for (i = 0; i < im->w; i++)
  126.             pp1[j][i] = getc(fp);
  127.         }
  128.       }
  129.     else
  130.       {
  131.       for (j = 0; j < im->h; j++)
  132.         {
  133.         REPORT(j, rlines);
  134.         for (i = 0; i < im->w; i++)
  135.           {
  136.               err += ((fcc = readpint(fp)) < 0);
  137.               if (fcc > cmax)
  138.               fcc = cmax;
  139.               if (fcc < cmin)
  140.               fcc = cmin;
  141.               if (cmax - PCMAXV)
  142.               fcc *= ((float) PCMAXV / cmax);
  143.               pp1[j][i] = fcc;
  144.           }
  145.         if (err)
  146.           {
  147.               fprintf(stderr, "total %d errors encountered\n", err);
  148.               remove_progress_report();
  149.               return -1;
  150.           }
  151.         }
  152.       }
  153.     im->mraster = pp1;
  154.  
  155.     get_bit_defmap(lindex, im->cmap->ct[0],
  156.            im->cmap->ct[1],
  157.            im->cmap->ct[2],
  158.            im->cmap->colors);
  159.     im->type = lindex ? T_CMAP : T_GMAP;
  160.     (void) fill_image_struct(im, pp1, im->h, im->w, im->type);
  161.     /* misc. info. All starts with @ */
  162.     xyequal = yx;
  163.     mag = 0;
  164.     t2b = 1;
  165.     while (fgets(tmpstr, sizeof(tmpstr), fp))
  166.       {
  167.       if (strncmp(tmpstr, "@MAP", 4) == 0)
  168.         {
  169.         int nentries;
  170.         err = 0;
  171.         sscanf(tmpstr + 5, "%d", &nentries);
  172.         holo_getmap(im, fp, nentries);
  173.         }
  174.       else if (strncmp(tmpstr, "@ASP", 4) == 0)
  175.         {
  176.         xyequal = 1;
  177.         }
  178.       else if (strncmp(tmpstr, "@YX", 3) == 0)
  179.         {
  180.         yx = 1;
  181.         }
  182.       else if (strncmp(tmpstr, "@B2T", 4) == 0)
  183.         {
  184.         t2b = 0;
  185.         }
  186.       else if (strncmp(tmpstr, "@MAG", 4) == 0)
  187.         {
  188.         mag = 1;
  189.         }
  190.       else if (strncmp(tmpstr, "@XUNIT", 6) == 0)
  191.         {
  192.         xunit = atof(tmpstr + 7);
  193.         }
  194.       else if (strncmp(tmpstr, "@YUNIT", 6) == 0)
  195.         {
  196.         yunit = atof(tmpstr + 7);
  197.         }
  198.       else if (strncmp(tmpstr, "@X0", 3) == 0)
  199.         {
  200.         xori = atoi(tmpstr + 4);
  201.         }
  202.       else if (strncmp(tmpstr, "@Y0", 3) == 0)
  203.         {
  204.         yori = atoi(tmpstr + 4);
  205.         }
  206.       else if (strncmp(tmpstr, "#T", 2) == 0)
  207.         {            /* text */
  208.         fseek(fp, -(long) strlen(tmpstr), SEEK_CUR);
  209.         load_text(im);
  210.         }
  211.       else if (strncmp(tmpstr, "#S", 2) == 0)
  212.         {            /* text */
  213.         fseek(fp, -(long) strlen(tmpstr), SEEK_CUR);
  214.         load_sgf(im);
  215.         }
  216.       }
  217.  
  218.     set_marking_units(1.0 / xunit, 1.0 / yunit);
  219.     set_marking_origin(xori, yori);
  220.  
  221.     /* we are officially done. Now process options */
  222.     pp = 0;
  223.     if (mag)
  224.       {                /* always keep aspect ratio */
  225.       xfac = def_w / im->w;
  226.       yfac = def_h / im->h;
  227.       factor = (xfac > yfac) ? yfac : xfac;
  228.       nw = im->w * (factor ? factor : 1);
  229.       nh = im->h * (factor ? factor : 1);
  230.       if (nw - im->w || nh - im->h)
  231.           img_scale(im, nh, nw, 0, 0, 0);
  232.       }
  233.     else
  234.       {
  235.       if (xyequal || im->w == im->h)
  236.         {
  237.         nw = (im->w < w_min) ? def_w : im->w;
  238.         nh = (im->h < h_min) ? def_h : im->h;
  239.         }
  240.       else
  241.         {
  242.         xfac = def_w / im->w;
  243.         yfac = def_h / im->h;
  244.         factor = (xfac > yfac) ? yfac : xfac;
  245.         nw = im->w * (factor ? factor : 1);
  246.         nh = im->h * (factor ? factor : 1);
  247.         }
  248.       if (nw - im->w || nh - im->h)
  249.         {
  250.         /*
  251.          * could've been using scale_hifi, but we want the
  252.          * interpolation to happen for the index ....
  253.          */
  254.         pp = interpol2(pp1, im->h, im->w, nh, nw, 1);
  255.         }
  256.       }
  257.     remove_progress_report();
  258.     if (pp || pp1 != im->mraster)
  259.       {
  260.       sprintf(tmpstr, "(was %-dX%-d)", im->w, im->h);
  261.       strcpy(im->misc, tmpstr);
  262.       }
  263.     return pp ? fill_image_struct(im, pp, nh, nw, im->type) : 0;
  264. }
  265.  
  266. /* ARGSUSED */
  267. const char *
  268. HOLO_wdefault(const IPTR im)
  269. {
  270.     return raw ? "Raw" : "Ascii";
  271. }
  272.  
  273. /* ARGSUSED */
  274. int
  275. HOLOdump_init(IPTR im)
  276. {
  277.     raw = !raw;
  278.     return 0;
  279. }
  280.  
  281. /* write to disk */
  282. int
  283. HOLO_dump(IPTR im)
  284. {
  285.     register ci_t *p = im->raster, *q;
  286.     register int i;
  287.     register pc_t *gr = im->cmap->ct[3];
  288.     int gmap = 0;        /* always zero for the moment BUG */
  289.     FILE *fp = im->fp;
  290.     int ckk = (PCMAX < 1000);
  291.     int per_line = ckk ? 16 : 14;
  292.     const char *fmt = ckk ? "%4d" : "%5d";
  293.  
  294.     fprintf(fp, "%s\n%d %d\n 1 %d\n 1 %d\n 0 %d\n",
  295.         raw ? "DSBB" : "DSAA", im->w, im->h, im->w, im->h, (PCMAX - 1));
  296.  
  297.     show_busy("Writing ...");
  298.     pack_cmap_g(im->cmap);
  299.     if (!raw)
  300.       {
  301.       for (i = 0, q = p + im->w * im->h; p < q; p++, i++)
  302.         {
  303.         fprintf(fp, fmt, gmap ? (int) gr[*p] : (int) *p);
  304.         if ((i + 1) % per_line == 0)
  305.             putc('\n', fp);
  306.         }
  307.       if (i % per_line != 0)
  308.           putc('\n', fp);
  309.       }
  310.     else
  311.       {
  312.       for (q = p + im->w * im->h; p < q; p++)
  313.         {
  314.         putc(gmap ? gr[*p] : *p, fp);
  315.         }
  316.       }
  317.     end_busy();
  318.     return fflush(fp);
  319. }
  320.  
  321. static void
  322. holo_getmap(IPTR im, FILE * fp, int n)
  323. {
  324.     int i, j, c, err = 0;
  325.     unsigned char cb[PCMAX];
  326.  
  327.     if (n > PCMAXV || n <= 0)
  328.       {
  329.       M_err("HoloMap", "Bad number of colors");
  330.       return;
  331.       }
  332.     err = 0;
  333.     if (!raw)
  334.       {
  335.       for (j = 0; j < 3; j++)
  336.         {
  337.         for (i = 0; i < n; i++)
  338.           {
  339.               err += (fscanf(fp, "%d", &c) != 1);
  340.               im->cmap->ct[j][i] = c;
  341.           }
  342.         }
  343.       }
  344.     else
  345.       {
  346.       for (j = 0; j < 3; j++)
  347.         {
  348.         err += (fread(cb, 1, n, fp) != n);
  349.         for (i = 0; i < n; i++)
  350.             im->cmap->ct[j][i] = cb[i];
  351.         }
  352.       }
  353.     if (err)
  354.       {
  355.       M_err("HoloMap", "Error readin colormaps\n", stderr);
  356.       return;
  357.       }
  358.     return;
  359. }
  360.  
  361. /* cycle thru local maps */
  362. void
  363. HOLO_special(IPTR im)
  364. {
  365.     if (IS_CI(im))
  366.       {
  367.       ++lindex;
  368.       lindex %= get_number_of_defmaps();
  369.       get_bit_defmap(lindex, im->cmap->ct[0],
  370.              im->cmap->ct[1], im->cmap->ct[2], im->cmap->colors);
  371.       im->type = lindex ? T_CMAP : T_GMAP;
  372.       set_cmap(im->cmap);
  373.       }
  374. }
  375.  
  376.  
  377. /*
  378.  * general 2-D linear interpolation where old and new grids are mismatched,
  379.  * i.e., (or+1)/in != integer
  380.  */
  381. void *
  382. interpol2(ci_t **in, int ir, int ic, int or, int oc, int rep)
  383. {
  384.     register int i, j, kx, ky, kx1, ky1;
  385.     register float dx, dy, f1, f2, f, sx, sy;
  386.     register float *flookup;
  387.     register float tmp;
  388.     register ci_t **m;
  389.     int rlines;
  390.  
  391.     if (ir <= 1 || ic <= 1 || or <= 1 || oc <= 1 || !in)
  392.       {
  393.       fputs("Bad parameter in interpol\n", stderr);
  394.       return 0;
  395.       }
  396.     if ((m = get_mat(or, oc, sizeof(ci_t))) == 0)
  397.       return 0;
  398.     if ((flookup = malloc(sizeof(float) * (oc + 1))) == 0)
  399.       return 0;
  400.     sx = (float) (ir - 1) / (or - 1);
  401.     sy = (float) (ic - 1) / (oc - 1);
  402.     for (i = 0; i < oc; i++)
  403.     flookup[i] = sy * i;
  404.  
  405.     rlines = rep ? progress_report("Interpolating ...", or) : 10000;
  406.     /* f(x+dx,y+dy)=f(x,y)+f'|y dx + f'|x dy + f'' dxdy */
  407.     for (i = 0; i < or; i++)
  408.       {
  409.       if (rep)
  410.         {
  411.         REPORT(i, rlines);
  412.         }
  413.       dx = sx * i;
  414.       kx = dx;
  415.       kx1 = (kx == ir - 1) ? kx : kx + 1;
  416.       dx = dx - kx;
  417.       for (j = 0; j < oc; j++)
  418.         {
  419.         ky = flookup[j];
  420.         ky1 = ky + (ky < ic - 1);
  421.         dy = flookup[j] - ky;
  422.         tmp = (float) in[kx][ky];
  423.         f1 = (float) in[kx1][ky] - tmp;
  424.         f2 = (float) in[kx][ky1] - tmp;
  425.  
  426.         f = tmp + dx * f1 + dy *
  427.             (f2 + dx * (in[kx1][ky1] - tmp - f1 - f2)) + 0.6;
  428.         m[i][j] = (ci_t) f;
  429.         }
  430.       }
  431.     free(flookup);
  432.     remove_progress_report();
  433.     return m;
  434. }
  435.  
  436. #endif
  437.